The present invention relates to an image processing method, an image processing apparatus, an image processing program and an image recording apparatus; particularly to an image processing method, an image processing apparatus, an image processing program and an image recording apparatus for an image formed into image signals through scanning of a color photographic film.
The image formed on the color photographic film is optically read by a CCD sensor or the like and is converted into image signals. After having been subjected to various types of image processing represented by negative/positive reversal, brightness adjustment, color balance adjustment, removal of granular noise and enhancement of sharpness, such image signals are distributed through such media as a CD-R, floppy (R) disk and memory card or via the Internet, and are outputted as hard copy images on silver halide photographic paper by an inkjet printer, thermal printer or the like. Alternatively, such image signals are displayed on the medium such as CRT, liquid crystal display or plasma display to be viewed.
Generally, a photographic film image is formed by a collection of pigment clouds of various sizes. When this image is enlarged for observation, mottled granular irregularity is visible according to the size of pigment clouds, although uniform colors should appear. Corresponding granular noise signals are included in the image signals obtained by optical reading of the image formed on a photographic film, using a CCD sensor or others. These granular noise signals are increased considerably with image processing of sharpness enhancement in particular, with the result that image quality is deteriorated. This has created a problem.
In recent years, a less costly digital still camera (hereinafter abbreviated as “DSC”) has come into widespread use. The DSC incorporated in such equipment as a cellular phone and laptop PC is also extensively used. The image sensor used in a less-costly DSC is characterized by a small pixel pitch. Shot noise tends to be produced at a low sensitivity, and not much consideration is given to cooling of an image sensor, so that conspicuous dark current noise is produced. The CMOS image sensor is often adopted in the less-costly DSC, so leakage current noise is conspicuous. When such noise is further subjected to image processing of interpolation of color filter arrangement and edge enhancement, the mottled granular irregularities are formed to deteriorate image quality. This has raised a problem. (For DSC noise and interpolation of color film arrangement, see the Non-Patent Document 1 listed later, for example).
To remove noise from image signals, the low-pass filter median filter technique is known. (See the Non-Patent Document 2 listed later, for example). However, noise removal by simple filtering involves reduced image sharpness, and a satisfactory image cannot be obtained.
Another widely known method for solving the above-mentioned problem is to use multiple filters to separate image signals into multiple frequency band components, and then to carry out suppression and enhancement for each frequency band. The Patent Documents 1 and 2 (listed later) propose the following technique, for example: Inputted image signals are decomposed into low, intermediate and high frequency band components, and enhancement processing is applied to the above-mentioned high frequency band components. At the same time, suppression processing is applied to the above-mentioned intermediate frequency band components. Processed high and intermediate frequency band components are synthesized with low frequency band components to get processed image signals, whereby granular noise is suppressed and sharpness is enhanced. By uniform suppression of the intermediate frequency band component where granular noise is mainly present, this technique has the effect of removing the granular irregularities that appear as mottles on the screen. However, since the information on image structure present on the intermediate frequency band components is also suppressed at the same time, shadow is suppressed on the bridge of the nose or around the eyes, for example, and the face appears blurred. A smooth expressionless face wearing makeup will appear.
The Patent Documents 3 and 4 (listed later) propose a technique of changing the filter conditions for each local site of an image, using the low-pass filter with varying sizes and shape. In this technique, the intermediate and high frequency band components are suppressed to remove the granular noise on the flat portion characterized by monotonous changes in brightness such as the cheek of the face or blue sky; whereas, on the edge portion such as hair and contour of the face characterized by sharp change in the brightness, a blur can be avoided without the intermediate frequency band components being suppressed. To get such an effect, however, a decision must be made to set the size and shape of the low-pass filter adequately for each step of processing. This takes much time and it is difficult to get sufficient effects by automatic processing of a desired image.
The frequency band is divided for each local site of an image, and a technique of using a wavelet transform is known as an effective way of suppression and enhancement for each frequency band. The details of wavelet transform are disclosed in the Non-Patent Documents 3 and 4. The following describes the overview:
The wavelet transform is operated as follows: In the first place, the following wavelet function is used, where vibration is observed in a finite range as shown in FIG. 1:
                                                        [                              Eq                .                                                                  ⁢                1                            ]                                                                                                            ψ                                      a                    ,                    b                                                  ⁡                                  (                  x                  )                                            =                              ψ                ⁢                                                                  ⁢                                  (                                                            x                      -                      b                                        a                                    )                                                                                        (        1        )            
Using the above function, the wavelet transform coefficient <f, ψa, b>with respect to input signal f(x) is obtained by:
                                                        [                              Eq                .                                                                  ⁢                2                            ]                                                                                          〈                                  f                  ,                                      ψ                                          a                      ,                      b                                                                      〉                            ≡                                                1                  a                                ⁢                                  ∫                                                                                    f                        ⁡                                                  (                          x                          )                                                                    ·                                              ψ                        ⁡                                                  (                                                                                    x                              -                              b                                                        a                                                    )                                                                                      ⁢                                          ⅆ                      x                                                                                                                              (        2        )            
Through this process, input signal is converted into the sum total of the wavelet function.
                                                        [                              Eq                .                                                                  ⁢                3                            ]                                                                                          f                ⁡                                  (                  x                  )                                            =                                                ∑                                      a                    ,                    b                                                  ⁢                                                      〈                                          f                      ,                                              ψ                                                  a                          ,                          b                                                                                      〉                                    ·                                                            ψ                                              a                        ,                        b                                                              ⁡                                          (                      x                      )                                                                                                                              (        3        )            
In the above equation, “a” denotes the scale of the wavelet function, and “b” the position of the wavelet function. As shown in FIG. 1, as the value “a” is greater, the frequency of the wavelet function ψa, b(x) is smaller. The position where the wavelet function ψa, b(x) vibrates moves according to the value of position “b”. Thus, Eq. 3 signifies that the input signal f(x) is decomposed into the sum total of the wavelet function ψa, b(x) having various scales and positions.
A great number of the wavelet functions are known, that allow the above-mentioned conversion. In the field of image processing, orthogonal wavelet and biorthogonal wavelet biorthogonal wavelet are put into common use. The following describes the overview of the conversion calculation of the orthogonal wavelet and biorthogonal wavelet.
Orthogonal wavelet and biorthogonal wavelet functions are defined as follows:
                                                        [                              Eq                .                                                                  ⁢                3                            ]                                                                                                            ψ                                      i                    ,                    j                                                  ⁡                                  (                  x                  )                                            =                                                2                                      -                    i                                                  ⁢                                  ψ                  ⁡                                      (                                                                  x                        -                                                  j                          ·                                                      2                                                                                                                                                      ⁢                              i                                                                                                                                                                            2                          ⁢                                                                                                                                i                                                              )                                                                                                          (        4        )            
where “i” denotes a natural number.
Comparison between Eq. 4 and Eq. 1 shows that the value of scale “a” is defined discretely by an i-th power of “2”, according to orthogonal wavelet and biorthogonal wavelet. This value “i” is called a level. In practical terms, level “i” is restricted up to finite upper limit N, and input signal is converted as follows:
                                                        [                              Eq                .                                                                  ⁢                5                            ]                                                                                                                                                                              f                        ⁡                                                  (                          x                          )                                                                    ≡                                            ⁢                                              S                        0                                                              =                                                                                            ∑                          j                                                ⁢                                                                              〈                                                                                          S                                0                                                            ,                                                              ψ                                                                  1                                  ,                                  j                                                                                                                      〉                                                    ·                                                                                    ψ                                                              1                                ,                                j                                                                                      ⁡                                                          (                              x                              )                                                                                                                          +                                                                        ∑                          j                                                ⁢                                                                              〈                                                                                          S                                0                                                            ,                                                              ψ                                                                  1                                  ,                                  j                                                                                                                      〉                                                    ·                                                                                    ϕ                                                              1                                ,                                j                                                                                      ⁡                                                          (                              x                              )                                                                                                                                                                                                                                            ≡                                        ⁢                                                                                            ∑                          j                                                ⁢                                                                                                            W                              1                                                        ⁡                                                          (                              j                              )                                                                                ·                                                                                    ψ                                                              1                                ,                                j                                                                                      ⁡                                                          (                              x                              )                                                                                                                          +                                                                        ∑                          j                                                ⁢                                                                                                            S                              1                                                        ⁡                                                          (                              j                              )                                                                                ·                                                                                    ϕ                                                              1                                ,                                j                                                                                      ⁡                                                          (                              x                              )                                                                                                                                                                                                                                  (        5        )                                                                                    S                                  i                  -                  1                                            =                            ⁢                                                                    ∑                    j                                    ⁢                                                            〈                                                                        S                                                      i                            -                            1                                                                          ,                                                  ψ                                                      1                            ,                            j                                                                                              〉                                        ·                                                                  ψ                                                  1                          ,                          j                                                                    ⁡                                              (                        x                        )                                                                                            +                                                      ∑                    j                                    ⁢                                                            〈                                                                        S                                                      i                            -                            1                                                                          ,                                                  ψ                                                      1                            ,                            j                                                                                              〉                                        ·                                                                  ϕ                                                  1                          ,                          j                                                                    ⁡                                              (                        x                        )                                                                                                                                                                    ≡                            ⁢                                                                    ∑                    j                                    ⁢                                                                                    W                        1                                            ⁡                                              (                        j                        )                                                              ·                                                                  ψ                                                  1                          ,                          j                                                                    ⁡                                              (                        x                        )                                                                                            +                                                      ∑                    j                                    ⁢                                                                                    S                        1                                            ⁡                                              (                        j                        )                                                              ·                                                                  ϕ                                                  1                          ,                          j                                                                    ⁡                                              (                        x                        )                                                                                                                                                    (        6        )                                                      f            ⁡                          (              x              )                                ≡                    ⁢                      S            0                          =                                            ∑                              i                =                1                            N                        ⁢                                          ∑                j                            ⁢                                                                    W                    i                                    ⁡                                      (                    j                    )                                                  ·                                                      ψ                                          1                      ,                      j                                                        ⁡                                      (                    x                    )                                                                                +                                    ∑              j                        ⁢                                                            S                  N                                ⁡                                  (                  j                  )                                            ·                                                ϕ                                      1                    ,                    j                                                  ⁡                                  (                  x                  )                                                                                        (        7        )            
The second term of Ex. 5 denotes that the low frequency band component of the residue that cannot be represented by the sum total of wavelet function ψ1, j(x) of level 1 is represented in terms of the sum total of scaling function φ1, j(x). An adequate scaling function in response to the wavelet function is employed (See Non-patent Documents 3 and 4 listed later). This means that input signal f(x)≡S0 is decomposed into the high frequency band component W1 and low frequency band component Si of level 1 by the wavelet transform of level 1 shown in Eq. 5. Since the wavelet function ψi, j(x) of the minimum traveling unit of the wavelet function ψi, j(x) is 2i, each of the signal volume of high frequency band component W1 and low frequency band component S1 with respect to the signal volume of input signal “S0” is ½. The sum total of the signal volumes W1 and S1 is equal to the signal volume of input signal “S0”. The low frequency band component S1 of level 1 is decomposed into high frequency band component W2 and low frequency band component S2 of level 2 by Eq. 6. After that, transform is repeated up to level N, whereby input signal “S0” is decomposed into the sum total of the high frequency band components of levels 1 through N and the sum of the low frequency band components of level N, as shown in FIG. 7.
Here the wavelet transform of level 1 shown in Eq. 6 is known to be computed by filtering, as shown in FIG. 2 (See Non-Patent Documents 3 and 4). In FIG. 2, LPF denotes a low-pass filter and HPF a high-pass filter. An appropriate filter coefficient is determined in response to the wavelet function (See Non-Patent Documents 3 and 4). Symbol 2↓ shows the down sampling where every other samples are removed (thinned out). The wavelet transform of level 1 in the secondary signal such as image signal is computed by the processing of filtering as shown in FIG. 3. In FIG. 3, LPFx, HPFx and 2↓x denote processing in the direction of “x”, whereas LPFy, HPFy and 2↓y denote processing in the direction of “y”. The low frequency band component Sn−1 is decomposed into three high frequency band components Wvn, Whn, Wdn and one low frequency band component Sn by the wavelet transform of level 1. Each of the signal volumes of Wvn, Whn, Wdn and Sn generated by decomposition is ½ that of the Sn−1 prior to decomposition in both vertical and horizontal directions. The total sum of signal volumes of four components subsequent to decomposition is equal to the signal Sn−1 prior to decomposition. FIG. 4 is a schematic diagram representing the process of the Input signal S0 being decomposed by the wavelet transform of level 3.
Further, when wavelet inverse transform is applied to Wvn, Whn, Wdn and Sn generated by decomposition, the signal Sn−1 prior to decomposition is known to be re-configured completely. In FIG. 5, LPF′ denotes a low-pass filter and HPF′ a high-pass filter. In the case of orthogonal wavelet, the same coefficient as that used in the wavelet transform is used as this filter coefficient; whereas in the case of biorthogonal wavelet, the coefficient different from that used in the wavelet transform is used as this filter coefficient. (See the above-mentioned Reference Documents). Further, 2↑ denotes the up-sampling where zero is inserted into every other signals. The LPF′x, HPF′x and 2↑x denote processing in the direction of “x”, whereas LPF′y, HPF′y and 2↓y denote processing in the direction of “y”.
The following image processing method is proposed as a known technique: Image signals are decomposed into multiple frequency band components by such a orthogonal wavelet and a biorthogonal wavelet, and each frequency band components is edited (data-processed). After that, inverse-transform is applied to configure image signals with reduced noise. According to Patent Document 5, the image signal representing the radiation screen is subjected to wavelet transform, whereby the above-mentioned image signals are decomposed into multiple frequency band components. The specified frequency band component out of the above-mentioned multiple frequency band components is subjected to the processing wherein signal values where the absolute value of each of the above-mentioned signal value is below the specified threshold value is reduced to 0. When inverse wavelet transform is applied to the frequency band component subjected to the above-mentioned processing and other frequency band components, image signals with reduced noise are obtained. In this image processing method, however, when the above-mentioned is applied to color images, the RGB balance close to the edge of the subject is lost, and a false color contour appears. This is very unseemly. In the noise structure of image signals gained by optical reading of the image formed on the color photographic film with a CCD sensor, mottled granular irregularities based on the size of the coloring pigment cloud is predominant, unlike the radiation image. If such a big threshold value as to erase this mottled irregularities is set, the sharpness of the image will be lost, or an artifice will appear. Conversely, when the threshold value is S0 small that the sharpness of the image is maintained, the mottled irregularities will not be erased. For these reasons, application of the technology of Patent Document 5 to a color image does not ensure satisfactory effects.
According to Tokkai 2000-224421, wavelet transform is applied to image signals to decompose the above-mentioned image signals into multiple frequency band components. Noise extraction processing is applied to the specified frequency band components, and noise elimination processing is applied to the above-mentioned specified frequency band components, based on the result of above-mentioned extraction processing, whereby processed frequency band components are obtained. Wavelet transform processing is applied to the above-mentioned processed frequency band components to get the low frequency band component one-step lower (of higher level) than the above-mentioned specified frequency band components. Processed frequency band component is obtained for each frequency band component by repeating the above-mentioned extraction processing wherein the above-mentioned low frequency band component one-step lower is the above-mentioned specified frequency band component, the above-mentioned noise elimination processing and the above-mentioned wavelet transform processing, until a desired frequency band is reached. Image signals with reduced noise are obtained by applying inverse wavelet transform to the above-mentioned processed frequency band component. However, when the above-mentioned technology is applied to color images, RGB balance is lost in the area where mottled irregularities are erased, with the result that false color spots appear. This is very unseemly. Further, noise extraction of low frequency signals and noise elimination processing must be repeated for each one-level transform, imposing a heavy computational load. Further, the medical image assumed in Patent Document 6 is restricted in the type of a subject, and comparatively monotonous images often occur. By contrast, the color photographic image is characterized by the mixed areas of different picture quality—an area where fine structures are closely packed in an image, the flat area, a bright area and a dark area, as in a people portrait against the background of a forest. In such a color image, noise elimination conditions must be changed, based on the area structure captured in perspective. According to the method given in the Patent Document 6, the noise elimination processing condition on the resolution level “n” must be determined according to the information on the above-mentioned resolution of level “n”. This makes it difficult to design an algorithm for determining the noise elimination conditions.
Non-Patent Documents 5, 6 and 7 give a detailed explanation to the Dyadic Wavelet used in the present invention. The following gives the overview:
The wavelet function of the Dyadic Wavelet is defined as follows:
                                                        [                              Eq                .                                                                  ⁢                6                            ]                                                                                                            ψ                                      i                    ,                    j                                                  ⁡                                  (                  x                  )                                            =                                                2                                      -                    i                                                  ⁢                                  ψ                  ⁡                                      (                                                                  x                        -                        j                                                                    2                                                                                                                                  ⁢                          i                                                                                      )                                                                                                          (        8        )            
where “i” denotes a natural number.
Wavelet functions of orthogonal wavelet and biorthogonal wavelet are discretely defined when the minimum traveling unit of the position on level “i” is 2i, as described above. By contrast, in the two-term wavelet, the minimum traveling unit of the position is constant, despite level “i”. This difference provides the Dyadic Wavelet transform with the following characteristics:
Characteristic 1: The signal volume of each of high frequency band component Wi and low frequency band component Si generated by the Dyadic Wavelet transform is the same as that of signal Si−1 prior to transform.
                                                        [                              Eq                .                                                                  ⁢                7                            ]                                                                                                                                                      S                                              i                        -                        1                                                              =                                        ⁢                                                                                            ∑                          j                                                ⁢                                                                              〈                                                                                          S                                                                  i                                  -                                  1                                                                                            ,                                                              ψ                                                                  i                                  ,                                  j                                                                                                                      〉                                                    ·                                                                                    ψ                                                              i                                ,                                j                                                                                      ⁡                                                          (                              x                              )                                                                                                                          +                                                                        ∑                          j                                                ⁢                                                                              〈                                                                                          S                                                                  i                                  -                                  1                                                                                            ,                                                              ϕ                                                                  i                                  ,                                  j                                                                                                                      〉                                                    ·                                                                                    ϕ                                                              i                                ,                                j                                                                                      ⁡                                                          (                              x                              )                                                                                                                                                                                                                                            ≡                                        ⁢                                                                                            ∑                          j                                                ⁢                                                                                                            W                              i                                                        ⁡                                                          (                              j                              )                                                                                ·                                                                                    ψ                                                              i                                ,                                j                                                                                      ⁡                                                          (                              x                              )                                                                                                                          +                                                                        ∑                          j                                                ⁢                                                                                                            S                              i                                                        ⁡                                                          (                              j                              )                                                                                ·                                                                                    ϕ                                                              i                                ,                                j                                                                                      ⁡                                                          (                              x                              )                                                                                                                                                                                                                                  (        9        )            
Characteristic 2: The following relationship is found between the scaling function φi, j(x) and wavelet function ψi, j(x):
                                                        [                              Eq                .                                                                  ⁢                8                            ]                                                                                                            ψ                                      i                    ,                    j                                                  ⁡                                  (                  x                  )                                            =                                                                    ∂                                                                                                                      ∂                    x                                                  ⁢                                                      ϕ                                          i                      ,                      j                                                        ⁡                                      (                    x                    )                                                                                                          (        10        )            
Thus, the high frequency band component Wi generated by the Dyadic Wavelet transform represents the first differential (gradient) of the low frequency band component Si.
Characteristic 3: With respect to Wi·γi (hereinafter referred to as “compensated high frequency band component) obtained by multiplying the coefficient γi (see the above-mentioned Reference Document on Dyadic Wavelet)) determined in response to the level “i” of the Wavelet transform, by high frequency band component, the relationship between levels of the signal intensities of compensated high frequency band components Wi·γi subsequent to the above-mentioned transform obeys a certain rule, in response to the singularity of the changes of input signals. To put it another way, the signal intensity of the compensated high frequency band component Wi·γi corresponding to smooth (differentiatable) signal changes shown by 1 and 4 of FIG. 6 increases with level number “i”; whereas the signal intensity of the compensated high frequency band component Wi·γi corresponding to stepwise signal changes shown by 2 of FIG. 6 stays constant independently of the level number “i”, and the signal intensity of the compensated high frequency band component Wi·γi corresponding to functional signal changes shown by 3 of FIG. 6 decreases with increase in level number “i”.
Characteristic 4: Unlike the above-mentioned method of orthogonal wavelet and biorthogonal wavelet, the method of Dyadic Wavelet transform on level 1 in the 2-D signals such as image signals is followed as shown in FIG. 7. The low frequency band component Sn−1 is decomposed into two high frequency band components Wxn, Wyn and one low frequency band component Sn by the wavelet transform of level 1. Two high frequency band components correspond to components x and y of the change vector Vn in the two dimensions of the low frequency band component Sn. The magnitude Mn of the change vector Vn and angle of deflection An are given by the following equation:
[Eq. 9]Mn=√{square root over (Wxn2+Wyn2)}  (11)An=argument (Wxn+iWyn)  (12)
It has been known that Sn−1 prior to transform can be re-configured when the Dyadic Wavelet inverse transform shown in FIG. 8 is applied to two high frequency band components Wxn, Wyn and one low frequency band component Sn.
In Patent Document 5, the following method is proposed to eliminate the white noise (Gaussian white noise) superimposed onto the monochrome image using the Dyadic Wavelet.
Step 1: The maximum of high frequency band component on each level is sought, and correspondence of the maximum positions between levels is established.
Step 2: When the absolute values of the associated maxima are reduced with increase in the level number, the maxima are removed.
Step 3: The linkage of remaining maxima on the plane surface is checked, and those having the linkage length in excess of the threshold value are preserved without being removed.
In addition to the above procedures, the following operations are also expounded:
Step 4: All signal values of the high frequency band component on the first level are discarded and signal values are synthesized according to the following method:                <1> The position of the maximum is made the same as that of the maximum on the second level left behind subsequent to the operations of steps 1 through 3.        <2> The magnitude M1 of the first level is determined by extrapolation of intensity between the levels higher than the second level of the magnitude Mn of the change vector defined in Eq. 11 in the above-mentioned position.        <3> The value A2 on the second level of the deflection angle An of the change vector defined in Eq. 12 in the above-mentioned position is copied to deflection angle A1 on the first level.        <4> The signal value of the high frequency band component on the first level is synthesized from the M1 and A1 obtained above.        
The method proposed above requires a great amount of computation, and much time and labor. Further, noise superimposed on image signals is assumed as white noise. S0 if this method is applied to the image containing mottled granular noise, as found in the case of the image gained by scanning a silver halide film or less costly DSC image, granular noise elimination will be insufficient or part of the granular noise will be enhanced. Further, if the above-mentioned method is applied to each of the RGB planes of a color image, a new fine granular noise similar to color misregistration will appear, raising another problem.
The documents cited in the above descriptions are listed as follow:
Patent Document 1: Tokkaihei 9-22460Patent Document 2: Tokkai 2000-215307Patent Document 3: Tokkai 2001-143068Patent Document 4: Tokkai 2001-155148Patent Document 5: Tokkaihei 9-212623Patent Document 6: Tokkai 2000-224421
Non-Patent Document 1: “Digital Photography” Chapter 2 and 3, published by The Society of Photographic Science and Technology of Japan, Corona Publishing Co., Ltd.
Non-Patent Document 2: “Practical Image Processing learnt in C-language” P54, by Inoue et al., Ohm Publishing Co., Ltd.
Non-Patent Document 3: “Wavelet and Filter Banks” by G. Strang & T. Nguyen, Wellesley-Cambridge Press
Non-Patent Document 4: “A wavelet tour of signal processing 2ed.” by S. Mallat, Academic Press
Non-Patent Document 5: “Singularity detection and processing with wavelets” by S. Mallat and W. L. Hwang, IEEE Trans. Inform. Theory 38 617 (1992)
Non-Patent Document 6: “Characterization of signal from multiscale edges” by S. Mallet and S. Zhong, IEEE Trans. Pattern Anal. Machine Intel. 14 710 (1992)
Non-Patent Document 7: “A wavelet tour of signal processing 2ed.” by S. Mallat, Academic Press